;This function creates a two dimensional komolgorov wavefront using fourier
;methods with random phases
;Size: the square array size
;L   : the outer scale in pixels (infinite if not set)
;seed: for the randomu algorithm (an unititialised variable 'seed' will do) 
;returns: a wavefront, which needs to be multiplied by some factor to get
;phase.
;The wavefront has an r_0 of 1 pix. To change this, divide by
;r_0^(5/6)
;if periodic, then the statistics will start to break down at the 10%
;level at size/50. if not periodic, then this will happen at size/25.

function kmf, sz, seed, L=L, periodic=periodic

size = sz
if (keyword_set(L) eq 0) then L = 1e6
if (keyword_set(periodic) eq 0) then periodic = 0
if periodic then size = size/2
ft_wf = complexarr(size*2,size*2)
phases = randomu(seed,size*2,size)
dist2 = dist(size*2, size*2)^2

;for i = 0,size-1 do begin
; for j = 0,size*2-1 do begin
;  if (i ne 0) or (j ne 0) then begin
;   dist2 = float(i)^2+float(((j+size) mod (size*2))-size)^2
;   ft_wf(i,j) = (dist2 + (size/(L*!pi))^2)^(-11.0/12.0)*exp(complex(0,6.283*phases(i,j)))
;  endif
;  ft_wf((size*2-i)mod (size*2),(size*2-j)mod (size*2)) = conj(ft_wf(i,j))
;  ;ft_wf(size*2-i,j) = ft_wf(i,j)
;  ;ft_wf(size*2-i,size*2-j) = ft_wf(i,j) 
; endfor
;endfor

;NB the 8.0 is a fudge factor... and the 0.1466 WAS 0.1513. There was
;no comment as to where this came from. The new number comes from 
;empirically checking a wavefront 4096x4096, with 10 pixel separation.
ft_wf[*,1:size-1] = 8.0*0.14667*((dist2[*,1:size-1] + (size/(L*!pi))^2)/float(size)^2)^(-11.0/12.0)*exp(complex(0,6.283*phases[*,1:size-1]))*float(size)
ft_wf[1:size-1,0] = 8.0*0.14667*((dist2[1:size-1,0] + (size/(L*!pi))^2)/float(size)^2)^(-11.0/12.0)*exp(complex(0,6.283*phases[1:size-1,0]))*float(size)

;for i = 0,size-1 do for j = 0,size*2-1 do $
;ft_wf((size*2-i)mod (size*2),(size*2-j)mod (size*2)) = conj(ft_wf(i,j))
ft_wf[1:*,size+1:2*size-1] = conj(reverse(reverse(ft_wf[1:*,1:size-1],1),2))
ft_wf[0,size+1:2*size-1] = conj(reverse( ft_wf[0,1:size-1],2 ))
ft_wf[size+1:2*size-1,0] = conj(reverse( ft_wf[1:size-1,0] ))
ft_wf[size, *] = complex(0.,0.) ;for symmetry. actually, if I could be bothered I'd
		   ;use the nyquist frequency too...

;mult = 2.41
if periodic then return, float(fft(ft_wf,-1));*mult
wf = extrac(fft(ft_wf,-1),0,0,size,size);*mult 

return, float(wf)
end
